1 Introduction

The purpose of this document is to provide a reproducible record of all analyses and figures in the main article. The main article is focused on the effect of response diversity on community stability in fluctuating environments. We are going to look at the effect of response diversity, richness, temperature and nutrients on community temporal stability. Specifically, we are going to look at the effect of fundamental balance (our measurement of stability) on temporal stability. Then, as response diversity is thought to stabilize temporal stability of aggregate community properties via asynchrony, we are going to look at the relationship between response diversity and asynchrony. Finally, as multiple evidence suggests that compensatory dynamics and temporal stability are determine by species interactions, we are going to analyse the effect of species interactions on stability to understand if they are more important than response diversity in driving temporal stability of community biomass.

This document is produced by an Rmarkdown file that includes code to reproduce from data all results presented in the main article.

2 Load datasets, Data wrangling and balance calculation

divergence_df <- read_csv("Data/divergence_df.csv")
load("Data/dens_biomass_poly.RData")

dd_all_pred<-read.csv("Data/morph_dd_pred.csv")
dd_all_pred_nonoise<-read.csv("Data/morph_dd_pred_nonoise.csv")

load("Data/ciliate_traits.Rdata")

df_slopes <- read_csv("Data/df_slopes_cor.csv")

# needs to have id_new variable
ciliate_traits <- ciliate_traits %>%
  dplyr::mutate(
    # Remove dots from the date
    cleaned_date = gsub("\\.", "", date),
    # Extract the part of id after the underscore
    id_suffix = sub(".*_(.*)", "\\1", id),
    # Combine cleaned_date, id_suffix, and species_initial into a new variable
    id_new = paste0(cleaned_date, id_suffix, composition)
  ) %>%
  # Optionally, remove the intermediate columns to clean up
  dplyr::select(-cleaned_date, -id_suffix,-new_id)

uniqueN(ciliate_traits$id_new)==nrow(ciliate_traits) # all unique  ;)

id_dd<-full_join(dd_all_pred,dplyr::select(ciliate_traits,id_new,biomass),join_by("id_new"))


## add day variable

#create a day variable from the date variable

id_dd<-dplyr::mutate(id_dd,date=as.Date(date,format = "%d.%m.%y"))

earliest_date<-min(id_dd$date)
days_since_earliest<-as.numeric(id_dd$date-earliest_date)+1
id_dd<-id_dd%>%dplyr::mutate(day=days_since_earliest)

#create a summarised df on microcosm level with each species seperate
# Make sure, that we have n_frames and not N_frames
names(id_dd)[names(id_dd) == "N_frames"] <- "n_frames"

#extrapolation_factor <- 9.301902  # for 16 x magnification 
extrapolation_factor <- 9.828125  # for 25 x magnification 
video_biomass_species <- c( "C", "P", "S","D","L","T")

biomasses <- id_dd %>%
  dplyr::group_by( day,temperature,nutrients,sample_ID,composition,predict_spec) %>% # group  by xxx
  dplyr::summarize(
    biomass = sum(biomass * n_frames, na.rm = TRUE) / (1 * 125) # if not 3 videos corrections is done below with dens_factor
  ) %>%
  dplyr::mutate(
    biomass = biomass * extrapolation_factor,
    )

biomasses<-biomasses%>%dplyr::mutate(biomass=biomass*1000)

dd_ts_id<-biomasses

#fill up missing dates with biomass<-0

fill_dd<-expand.grid(sample_ID=unique(dd_ts_id$sample_ID),day=unique(dd_ts_id$day),predict_spec=unique(dd_ts_id$predict_spec))
complete_ts<-full_join(fill_dd,dd_ts_id,join_by(sample_ID,day,predict_spec))

complete_ts$biomass[is.na(complete_ts$biomass)]<-0
complete_ts<-complete_ts%>%dplyr::mutate(composition=sub("_.*", "", sample_ID))
complete_ts<-complete_ts %>%
  dplyr::mutate(temperature = sapply(strsplit(as.character(sample_ID), "_"), function(x) paste(x[3], x[4], sep = "-")))
complete_ts<- dplyr::mutate(complete_ts,nutrients = gsub(".*Nut(.*?)_.*", "\\1", sample_ID))

# Now remove wrong combinations of composition and predict_spec / predict_spec

complete_ts<- complete_ts %>%
  rowwise() %>%
  dplyr::filter(predict_spec %in% unlist(strsplit(composition, ""))) %>%
  ungroup()  
complete_ts<-dplyr::mutate(complete_ts,temperature=as.character(temperature),
                    nutrients=as.character(nutrients),
                    richness=nchar(composition))

complete_ts<-complete_ts%>%group_by(sample_ID,composition,day)%>%dplyr::mutate(tot_biomass=sum(biomass))
complete_ts<-complete_ts%>%dplyr::mutate(biom_contribution=biomass/tot_biomass)

df_biomass_mod <- complete_ts

complete_ts<-complete_ts%>%dplyr::mutate(temperature=paste0(temperature," °C"),
                                      nutrients=paste0(nutrients," g/L"))


# introduce slopes of 
names(df_slopes)[names(df_slopes)=="species_initial"]<-"predict_spec"

slope_ts<-full_join(dplyr::select(df_slopes,nutrients,predict_spec,temperature,slope),complete_ts)
slope_ts<-slope_ts%>%dplyr::mutate(w_slope=biom_contribution*slope,
                            sign=sign(slope))

slope_ts<-slope_ts%>%group_by(sample_ID,temperature,nutrients,richness,composition,day,tot_biomass)%>%dplyr::summarize(
  sum_w_slopes=abs(sum(w_slope)),
                   mean_abs_slope=mean(abs(slope)),
  sum_abs_slope=sum(abs(slope)),
  abs_sum_slope=abs(sum(slope)),
  symmetry=abs(sum(sign)))


slope_ts<-slope_ts%>%dplyr::mutate(richness=as.factor(richness))


##create new variable where it checks, where the last observation =0 is; with complete_ts
aggr_ts <- slope_ts %>%
  group_by( sample_ID) %>%
  arrange(day) %>%
  mutate(
    # Create a flag for non-zero tot_biomass
    non_zero_biomass = tot_biomass != 0,
    # Find the last non-zero day
    last_non_zero_day = ifelse(any(non_zero_biomass), max(day[non_zero_biomass], na.rm = TRUE), NA),
    # Find the first zero day after the last non-zero day
    first_zero_day = ifelse(
      !is.na(last_non_zero_day),
      min(day[!non_zero_biomass & day > last_non_zero_day], na.rm = TRUE),
      NA
    ),
    # Flag for days after the first zero day
    is_after_first_zero_day = ifelse(!is.na(first_zero_day), day > first_zero_day, FALSE)
  ) %>%
  ungroup()

aggr_ts<-aggr_ts%>%mutate(rep_var=sub("_[^_]+$", "", sample_ID))

biomass_ts<-aggr_ts%>%group_by(day,temperature,nutrients,richness)%>%summarize(tot_biom=mean(tot_biomass),se_tot_biom=sd(tot_biomass)/sqrt(as.numeric(length(tot_biomass))))

3 Biomass

Let’s have a look at the biomass dynamics in the different environmental treatments.

3.0.1 tot biomass plot

Figure 1 : Community total biomass during the experiment in different environmental treatments. Different color represent richness levels.

4 Main Results

We now look at the main results of the experiment. We are going to look first at the effect of richness, temperature and nutrients on community temporal stability. Then, we are going to look at the relationship between divergence (original response diversity metric) and temporal stability. Finally, we are going to look at the relationship between response diversity and temporal stability.

In the whole analysis, we calculated the temporal stability of total community biomass as the inverse of the coefficient of variation (ICV) (i.e. \(\frac{\sigma}{\mu}\)).

4.0.1 Effect of T, N and R

Figure 2: Effects of richness (a), temperature (b), and nutrients (c) on community total biomass temporal stability.

We can see that richness does not have a clear effect on community temporal stability, while stability was higher at lower temperature, and nutrients increased community temporal stability.

4.0.2 Effect of Divergence

We look at the relationship between divergence (our original response diversity metric) and stability

Figure 3: Relationship between Divergence and temporal stability of total community biomass.

Divergence is positively related to temporal stability, suggesting that response diversity promotes stability. However, the relationship between divergence and stability becomes weaker as richness increases. We think that this is due to divergence considering only the responses of the 2 most “responding” species. Thus, when species richness increases, disregarding the responses of the other species in the community except the 2 responding the most makes the relationship between response diversity and stability weaker.

This is why, after running the experiment, we developed another metric to measure response diversity, which we called balance, and that is presented in the main text of the publication. Balance has several desirable features that makes it a more suitable metric than divergence: Independence of richness, higher predictive power, and accounts for the responses of all species in the community (as opposed to divergence that accounts for only the 2 most “responding” species).

Here, we provide extensive evidence of why balance is a better metric to measure response diversity than divergence, and thus justifying focusing the analysis around balance.

5 Comparing Divergence and Balance

5.1 Predictive power of Divergence and Balance

We first compare how well divergence and balance predict stability (predictive power).

5.1.1 Balance

# 

mod1 <- lm(data=complete_aggr,log10(stability)~log10(balance_f))

# Check model assumptions
#check_model(mod1)

5.1.2 Divergence

mod2 <- lm(data=complete_aggr,log10(stability)~(divergence))

# Check model assumptions
#check_model(mod2)

Table 1: Comparison of model performance of divergence and balance as predictors of stability. Model 1 has balance as predictor and model 2 has divergence as predictor.

model	AIC	AICc	BIC	R2	R2_adjusted	RMSE	Sigma
1	-89.27328	-89.17286	-78.79409	0.1917679	0.1884142	0.1510344	0.1516599
2	-55.71579	-55.61538	-45.23661	0.0720796	0.0682293	0.1618316	0.1625017

A model with Balance as predictor performs better than one with divergence as predictor, and it explains more of the variance in stability than divergence.

Moreover, from Figure 3, it looks like divergence declines in performance as richness increases. Let’s test this analytically. To do than we build a linear model having stability as response variable and either log10(balance) or divergence as predictor for each richness level. We then extract the R squared of the models and their standardised estimates. (standardized estimates were calculated centering divergence and balance using the function scale()).

# getting model estimates for each richness level
lm_divergence_richness_E <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
    results = map(model, broom::tidy)
  ) %>%
  unnest(results) %>% dplyr::filter(term=="scale(divergence)") 


# getting model R squared for each richness level

lm_divergence_richness_R <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
    results = map(model, broom::glance)
  ) %>%
  unnest(results)

# getting model estimatesf or each richness level
lm_balance_richness_E <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
    results = map(model, broom::tidy)
  ) %>%
  unnest(results) %>% dplyr::filter(term=="scale(log10(balance_f))") 



# getting model R squared for each richness level
lm_balance_richness_R <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
    results = map(model, broom::glance)
  ) %>%
  unnest(results)

Figure 4: Performance comparison of divergence vs balance. In (a), the R squared of linear models for divergence and balance are shown for each richness level. In (b), the estimates of the linear models for divergence and balance are shown for each richness level.

We can see that the R squared of divergence as predictor of stability becomes smaller as richness increases, while the R squared of balance as predictor of stability does not (actually increases slightly).

5.2 Comparing unique explanatory power of balance and divergence

Now we build a linear model were stability is modeled as a function of balance and divergence. Then, we compared the variance explained by the full model compared to a model containing either only balance or only divergence.

5.2.1 Full model - balance and divergence

lm_div_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)+divergence)

# Check model assumptions
# check_model(lm_div_balance)

5.2.2 model with only divergence

lm_div <- lm(data=complete_aggr,log10(stability)~divergence)

# Check model assumptions
# check_model(lm_div)

5.2.3 model with only balance

lm_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f))

# Check model assumptions
# check_model(lm_balance)

5.2.4 Comparision full model vs divergence only and balance only

Table 2: Comparison of model performance of divergence, balance and both as predictors of stability. Model 1 has both balance and divergence as predictors, model 2 has divergence as predictor, and model 3 has balance as predictor.

model	AIC	AICc	BIC	R2	R2_adjusted	RMSE	Sigma
1	-88.97683	-88.80876	-75.00458	0.1974141	0.1907259	0.1505060	0.1514437
1	-55.71579	-55.61538	-45.23661	0.0720796	0.0682293	0.1618316	0.1625017
2	-89.27328	-89.17286	-78.79409	0.1917679	0.1884142	0.1510344	0.1516599

5.2.5 Comparision full model vs balance only

Table 3: Anova table: a model with both balance and divergence as predictors is not significantly different from a model with only balance as predictor.

anova1 <- anova(lm_div_balance,  lm_balance)

# Convert to tidy format
anova_tidy1 <- broom::tidy(anova1)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy1 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	df.residual	rss	DF	Sum of Squares	F Statistic	p-value
log10(stability) ~ log10(balance_f) + divergence	240	5.504447	NA	NA	NA	NA
log10(stability) ~ log10(balance_f)	241	5.543171	-1	-0.03872444	1.688429	0.195

5.2.6 Comparision full model vs divergence only and divergence only

Table 4: Anova table: a model with both balance and divergence as predictors is significantly better from a model with only divergence as predictor.

anova2 <- anova(lm_div_balance,  lm_div)


anova_tidy2 <- broom::tidy(anova2)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy2 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	df.residual	rss	DF	Sum of Squares	F Statistic	p-value
log10(stability) ~ log10(balance_f) + divergence	240	5.504447	NA	NA	NA	NA
log10(stability) ~ divergence	241	6.364040	-1	-0.8595933	37.47922	0.000

Overall, balance explains more of the variance in stability than divergence, and there is virtually no difference between a model containing only balance and the full model.

5.3 Interaction divergence and richness

Richness had to be transformed to numeric and to be centered to avoid collinearity with divergence

lm_rich_div <- lm(data=complete_aggr,log10(stability)~divergence*scale(as.numeric(richness)))

# check model assumptions
# check_model(lm_rich_div)

Table 5: Type III anova table of the model with divergence and richness as predictors of stability.

anova3 <- car::Anova(lm_rich_div, type = "III")

anova_tidy3 <- broom::tidy(anova3)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy3 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	Sum of Squares	DF	F Statistic	p-value
(Intercept)	11.033652088	1	442.55571734	0.000
divergence	0.807044347	1	32.37025122	0.000
scale(as.numeric(richness))	0.001236238	1	0.04958505	0.824
divergence:scale(as.numeric(richness))	0.249582101	1	10.01064605	0.002
Residuals	5.958668583	239	NA	NA

Divergence significantly interact with richness, suggesting that the relationship between divergence and stability changes with richness. While an ideal metric of response diversity should be independent of richness.

We repeat the same model using balance instead of divergence.

lm_rich_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)*scale(as.numeric(richness)))

# check model assumptions
# check_model(lm_rich_balance)

Table 6: Type III anova table of the model with balance and richness as predictors of stability.

anova4 <- car::Anova(lm_rich_balance, type = "III")

anova_tidy4 <- broom::tidy(anova4)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy4 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	Sum of Squares	DF	F Statistic	p-value
(Intercept)	9.54928736	1	414.779284	0.000
log10(balance_f)	1.26534818	1	54.961191	0.000
scale(as.numeric(richness))	0.02471247	1	1.073401	0.301
log10(balance_f):scale(as.numeric(richness))	0.04049552	1	1.758948	0.186
Residuals	5.50239554	239	NA	NA

Balance does not significantly interact with richness, suggesting that the relationship between balance and stability is stable across richness levels.

5.4 Variable importance

Finally, we assess variable importance using the relative importance of predictors in the full model. We use the package vip (https://cran.r-project.org/web/packages/vip/vignettes/vip.html) to calculate the relative importance of predictors in the full model. The function vip::vip for multiple linear regression, or linear models (LMs), uses the absolute value of the -statistic as a measure of VI. Motivation for the use of the associated 𝑡-statistic is given in Bring (1994) [https://www.tandfonline.com/doi/abs/10.1080/00031305.1994.10476059].

vip::vip(lm_div_balance)

Figure 5: Variable importance in the model including both balance and divergence as predictors of stability.

We believe that the extensive evidence here provided justifies focusing the analysis around balance, and not divergence, as a metric of response diversity. We will thus only look at balance for the rest of the analysis.

6 Effect RD

We are now going to look at how response diversity (balance) affected temporal stability of total community biomass. We are going to look at the relationship between fundamental balance (so based only on species response surfaces measured in monoculture), an realised balance (measured accounting for species contribution to balance).

This is fundamentally testing our most important hypothesis.

Figure 6: Effects of fundamental and realised response diversity (measured as balance) on total community biomass temporal stability.

We can see that balance is always negatively related to temporal stability, which means that response diversity promotes stability across richness levels. Interestingly, we see that there is little difference between fundamental and realised balance. Yet, as the richness increases, the relationship between realised balance and stability becomes steeper compared to fundamental balance.

But is the difference in the slope of fundamental and realised balance significant? We can test this using a linear model with interaction between balance and type (factor: fundamental or realised balance).

6.1 Balance: realised vs fundamental

# compare if the slope of fundamental and realised balance is significantly different for each richness level
# Fit the linear model with interaction
model_F_R <- lm(log10(1/CV) ~ log10(balance) * type, data = main_r_dd)

Table 7: Type III anova table of the model with interaction between balance and type (fundamental vs realised) as predictors of stability.

Term	Sum of Squares	DF	F Statistic	p-value
ANOVA Table for Linear Model
(Intercept)	7.76166605	1	325.480366	0.000
log10(balance)	1.17164323	1	49.132089	0.000
type	0.09227462	1	3.869476	0.050
log10(balance):type	0.04940699	1	2.071850	0.151
Residuals	11.49415887	482	NA	NA

No significant difference in the slope of fundamental and realised balance was found.

6.1.1 Balance: realised vs fundamental by richness level

Table 7: Linear model results for the interaction between balance and type (fundamental vs realised) as predictors of stability for richness level.

Richness	Term	Estimate	Std_Error	T_Value	P_Value
Regression Results for log10(balance):type (Realised)
2	log10(balance):typerealised	−0.007	0.058	−0.116	0.908
3	log10(balance):typerealised	−0.030	0.029	−1.042	0.299
4	log10(balance):typerealised	−0.064	0.036	−1.785	0.076

Even within each richness level, the slope of fundamental and realised balance is never significantly different.

7 Linear models

7.1 Model: Fundamental balance

First we analyze the effect of fundamental balance, temperature, nutrients and richness on biomass temporal stability using a linear model. balance was modelled as continuous variables, while richness, temperature and nutrients were modelled as categorical variables. balance and stability were log-transformed to meet the assumptions of linear models.

lm_full<-lm(data=complete_aggr,log10(stability)~log10(balance_f)+(richness)+nutrients+temperature)

# check model assumptions
# check_model(lm_full)

Table 8: Linear model results for the effects of balance, richness, nutrients, and temperature on community stability. Estimates are presented with 95% confidence intervals and p-values.

Predictor	Linear Regression Results	95% CI¹	Linear Regression Results
Predictor	Estimate	95% CI¹	p-value
log10(balance_f)	-0.05	-0.08, -0.02	<0.001
richness
richness3 - richness2	-0.04	-0.09, 0.00	0.065
richness4 - richness2	-0.01	-0.06, 0.03	0.8
richness4 - richness3	0.03	-0.01, 0.07	0.3
nutrients
(0.35 g/L) - (0.01 g/L)	0.18	0.14, 0.22	<0.001
(0.75 g/L) - (0.01 g/L)	0.21	0.17, 0.26	<0.001
(0.75 g/L) - (0.35 g/L)	0.03	-0.01, 0.08	0.2
temperature
(22-25 °C) - (18-21 °C)	-0.08	-0.12, -0.03	<0.001
(25-28 °C) - (18-21 °C)	-0.10	-0.16, -0.04	<0.001
(25-28 °C) - (22-25 °C)	-0.02	-0.08, 0.04	0.7
¹ CI = Confidence Interval

A linear model was fitted to examine the effects of resource balance, richness, nutrients, and temperature on community stability (measured as log₁₀(stability)).

Among the predictors, log₁₀(balance) showed a significant negative effect on stability (Estimate = -0.05, SE = 0.016, p< 0.001). This suggests that as balance increases (more balance), stability tends to decrease.

Richness did not have a significant effect on stability within the conditions of this study.

Nutrient concentration also had a significant positive effect on stability, with estimates for 0.35 g/L (Estimate = 0.18, SE = 0.019, p < 0.001) and 0.75 g/L (Estimate = 0.21, SE = 0.019, p < 0.001) indicating increased stability with higher nutrient levels, when compared to the baseline (0.01 g/L).

Finally, temperature regimes showed a significant effect on stability. Both 22–25 °C (Estimate = -0.08, SE = 0.019, p < 0.001) and 25–28 °C (Estimate = -0.10, SE = 0.02, p < 0.001) significantly reduced stability when compared to the baseline (18–21 °C).

In summary, our findings show that temporal stability is significantly influenced by response diversity (balance), nutrient concentration, and temperature, with higher nutrient concentrations enhancing stability and higher temperatures reducing it. However, species richness was not a significant determinant of stability within the conditions of this study.

Table 9: Type II anova table of the model with balance, richness, nutrients, and temperature as predictors of stability.

Term	Sum of Squares	DF	F Statistic	p-value
Type III ANOVA Table for Linear Model
(Intercept)	2.11729127	1	152.932234	0.000
log10(balance_f)	0.15572514	1	11.248048	0.001
richness	0.07402585	2	2.673449	0.071
nutrients	1.98098126	2	71.543272	0.000
temperature	0.32820230	2	11.853048	0.000
Residuals	3.25348971	235	NA	NA

7.1.1 Interaction between temperature and nutrients

We may expect and interactive effect of the environmental variables on stability. We thus build a linear model with interaction between temperature and nutrients. However, there is high collinearity between temperature and nutrients, which may affect the model results.

lm_full_int<-lm(data=complete_aggr,log10(stability)~log10(balance_f)+(richness)+nutrients*temperature)

# check model assumptions
check_model(lm_full_int)

(#fig:model_check_int)model check 1.

So we transformed nutrients and temperature to numeric, and transformed temperature regimes in values = 1, 2, 3. Then, we centered the variables to avoid collinearity with the interaction term.

# transform nutrients and temperature to numeric. For this the units need to be removed, and temperature regimes should be transformed in values = 1, 2, 3
complete_aggr_2<- complete_aggr %>%
  # Remove the units from the 'nutrients' and 'temperature' columns
  mutate(
    nutrients = as.numeric(gsub(" g/L", "", nutrients)),  # Convert nutrients to numeric
    temperature = gsub(" °C", "", temperature)            # Remove the unit but keep as character
  ) %>%
  # Convert temperature ranges to numeric codes using case_when
  mutate(
    temperature = case_when(
      temperature == "18-21" ~ 1,
      temperature == "22-25" ~ 2,
      temperature == "25-28" ~ 3,
      TRUE ~ NA_real_         # Handle unexpected values with NA
    )
  )




# Fit the linear model with interaction
lm_full_int<-lm(data=complete_aggr_2,log10(stability)~log10(balance_f)+(richness)+scale(nutrients)*scale(temperature))

# check model assumptions
check_model(lm_full_int)

(#fig:model_check_int2)model check 1.

Table 10: Type III anova table of the model with interaction between temperature and nutrients as predictors of stability.

Term	Sum of Squares	DF	F Statistic	p-value
Type III ANOVA Table for Linear Model
(Intercept)	2.07662552	1	135.208125	0.000
log10(balance_f)	0.06339546	1	4.127649	0.043
richness	0.07254258	2	2.361607	0.096
scale(nutrients)	1.71980335	1	111.975599	0.000
scale(temperature)	0.32873431	1	21.403738	0.000
scale(nutrients):scale(temperature)	0.02595802	1	1.690115	0.195
Residuals	3.62466104	236	NA	NA

No interaction between nutrients and temperature was found to significantly affect stability. We thus retain the simplest model without interaction term.

8 Asynchrony

Response diversity (aka balance) has been suggested as a mechanism that promotes temporal stability of community biomass by promoting species asynchrony.

We thus calculated the asynchrony index suggested by Gross et al. (2014)[https://www.journals.uchicago.edu/doi/epdf/10.1086/673915] to calculate the effect of asynchrony on temporal stability and to see how reponse diversity relate to asynchrony. The index ranges between -1 and 1, with -1 indicating perfect asyncrony and 1 being perfectly synchronous, and 0 indicating random variation.

8.0.1 Plot stability vs. Asynchrony Gross

Figure 8: Relationship between temporal stability and asynchrony (Gross) divided by nutrient level.

The Pearson’s correlation between asynchrony and stability is significant (estimate = -0.23, p < 0.001).

cor.test((-1*async_aggr$synchrony_Gross),async_aggr$stability)

## 
##  Pearson's product-moment correlation
## 
## data:  (-1 * async_aggr$synchrony_Gross) and async_aggr$stability
## t = 3.7927, df = 239, p-value = 0.0001888
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1153693 0.3539711
## sample estimates:
##       cor 
## 0.2382622

8.0.2 Plot Asynchrony Gross vs fundamental balance

Figure 9: Relationship between asynchrony (Gross) and fundamental balance divided by nutrient level.

The Pearson’s correlation between asynchrony and balance is significant (estimate = 18, p = 0.003).

cor.test((-1*async_aggr$synchrony_Gross),(async_aggr$balance_f))

## 
##  Pearson's product-moment correlation
## 
## data:  (-1 * async_aggr$synchrony_Gross) and (async_aggr$balance_f)
## t = -2.9796, df = 239, p-value = 0.003184
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3082462 -0.0644258
## sample estimates:
##        cor 
## -0.1892515

9 Population stability

The relationship between community stability and the stability of the individual populations that make up the community is a key question in community ecology. Importantly, community stability can result from low population stability, if populations fluctuate asynchronously, or from high population stability, if populations do not fluctuate much. Synthesis of the literature suggests diversity can have a positive or negantive effect on population stability (Campbell et al 2010)[https://nsojournals.onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0706.2010.18768.x] and (Xu et al 2021)[https://onlinelibrary.wiley.com/doi/full/10.1111/ele.13777].

Theoretical work has suggested that community stability is a product of two quantities: the (a)synchrony of population fluctuations, and an average species-level population stability that is weighted by relative abundance (Thibaut & Connolly 2013)[https://onlinelibrary.wiley.com/doi/full/10.1111/ele.12019].

Critically, a balance value close to zero can result from high response diversity, but also from high population stability (population biomass does not change largely over time). We want to look now at whether our new metric of balance can capture these two stabilising mechanisms.

Thus, we first calculate species-level population stability weighted by relative abundance.

–>

10 SEM

Finally, we use a structural equation model (SEM) to explore how stability is influenced by asynchrony, population stability, balance and, nutrient levels. In order to develop a hypothesis regarding the influence of stability, we have drawn on existing literature. This has enabled us to posit that stability is influenced by two key factors: asynchrony and population stability. In turn, these are influenced by balance and, in our particular case, by nutrient levels.

## lavaan 0.6-19 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
## 
##                                                   Used       Total
##   Number of observations                           220         241
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                 1.511       1.433
##   Degrees of freedom                                 3           3
##   P-value (Chi-square)                           0.680       0.698
##   Scaling correction factor                                  1.055
##     Satorra-Bentler correction                                    
## 
## Model Test Baseline Model:
## 
##   Test statistic                               628.549     624.481
##   Degrees of freedom                                 9           9
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.007
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.007       1.008
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.008
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)                435.602     435.602
##   Loglikelihood unrestricted model (H1)             NA          NA
##                                                                   
##   Akaike (AIC)                                -847.204    -847.204
##   Bayesian (BIC)                              -806.481    -806.481
##   Sample-size adjusted Bayesian (SABIC)       -844.509    -844.509
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.087       0.082
##   P-value H_0: RMSEA <= 0.050                    0.825       0.847
##   P-value H_0: RMSEA >= 0.080                    0.067       0.055
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.087
##   P-value H_0: Robust RMSEA <= 0.050                         0.831
##   P-value H_0: Robust RMSEA >= 0.080                         0.067
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.006       0.006
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                      Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   stability ~                                                             
##     asynchrny_Grss      0.181    0.016   11.078    0.000    0.181    0.373
##     pop_stability       1.013    0.034   29.470    0.000    1.013    0.924
##   asynchrony_Gross ~                                                      
##     log_balance_f      -0.084    0.034   -2.493    0.013   -0.084   -0.176
##     nutrients          -0.199    0.025   -7.994    0.000   -0.199   -0.469
##   pop_stability ~                                                         
##     log_balance_f      -0.064    0.008   -7.748    0.000   -0.064   -0.301
##     nutrients           0.124    0.007   18.314    0.000    0.124    0.661
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .stability         0.217    0.017   12.489    0.000    0.217    1.333
##    .asynchrny_Grss   -0.110    0.060   -1.817    0.069   -0.110   -0.327
##    .pop_stability    -0.682    0.016  -41.719    0.000   -0.682   -4.596
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .stability         0.005    0.001    9.841    0.000    0.005    0.188
##    .asynchrny_Grss    0.088    0.012    7.616    0.000    0.088    0.781
##    .pop_stability     0.009    0.001    9.030    0.000    0.009    0.397
## 
## R-Square:
##                    Estimate
##     stability         0.812
##     asynchrny_Grss    0.219
##     pop_stability     0.603

Figure 10.1: SEM.

Model Fit Indices The model fit indices suggest that the model fits the data well.

Chi-Square Test (User Model): The chi-square test statistic for the user model is χ 2 =1.626 (scaled = 1.465) with 3 degrees of freedom and a p-value of 0.653 (scaled = 0.690). This indicates a good fit, as the test is non-significant, suggesting no significant difference between the observed and model-implied covariance matrices.

Comparative Fit Index (CFI) and Tucker-Lewis Index (TLI): Both CFI and TLI values are 1.000, indicating an excellent model fit. Values close to or above 0.95 are generally considered good.

Root Mean Square Error of Approximation (RMSEA): The RMSEA is 0.000, with a 90% confidence interval ranging from 0 to 0.090 (scaled = 0.080). This indicates a very good fit, as RMSEA values below 0.05 are ideal, and values below 0.08 are acceptable. The p-values for the RMSEA hypothesis tests suggest strong support for a close fit (RMSEA <= 0.05) and little evidence for a poor fit (RMSEA >= 0.08).

Standardized Root Mean Square Residual (SRMR): The SRMR value is 0.017, which is also within the acceptable range (values below 0.08 are generally considered good). Overall, the fit indices suggest that the model is an excellent fit for the data.

Regression Paths and Interpretation

Stability Regressions

Stability ~ Asynchrony_Gross (asynchrny_Grss): The standardized estimate for the effect of asynchrony on stability is 0.340 (p < 0.001), indicating a significant positive association. Higher asynchrony in species dynamics is associated with increased community stability.

Stability ~ Population Stability (pop_stability): The standardized estimate is 0.977 (p < 0.001), showing a strong positive relationship. This suggests that community stability is highly dependent on the stability of individual populations within the community.

Asynchrony_Gross Regressions

Asynchrony_Gross ~ Log10(Balance): The standardized estimate is -0.176 (p = 0.013), indicating a significant negative effect. Higher balance leads to lower asynchrony, suggesting that as balance increases, species within the community fluctuate more synchronously.

Asynchrony_Gross ~ Nutrients: The standardized estimate is -0.469 (p < 0.001), showing a strong negative relationship. Higher nutrient levels appear to reduce asynchrony, possibly by causing similar responses across species.

Population Stability Regressions

Population Stability ~ Log10(Balance): The standardized estimate is -0.296 (p < 0.001), indicating that higher balance is associated with lower population stability.

Population Stability ~ Nutrients: The standardized estimate is 0.635 (p < 0.001), showing that higher nutrient levels are associated with increased population stability, likely because nutrients enhance conditions that support stable population dynamics.

Variances and R-Squared Values R-Squared for Stability: The model explains 90.4% of the variance in community stability, indicating strong predictive power.

R-Squared for Asynchrony_Gross: The model explains 21.9% of the variance in asynchrony, which is moderate.

R-Squared for Population Stability: The model explains 56.2% of the variance in population stability, showing that nutrients and balance are important but not the only factors influencing it.

Summary Interpretation Model Fit: The model has an excellent fit, as indicated by the fit indices. Stability: Community stability is strongly influenced by both population stability and asynchrony among species, with population stability being the stronger predictor. Asynchrony and Balance: Asynchrony decreases with increasing balance and nutrients, suggesting that these factors promote more synchronized fluctuations among species. Population Stability and Nutrients: Higher nutrient levels are associated with increased population stability, suggesting that nutrient availability supports stable population dynamics. Conversely, higher balance is associated with decreased population stability.

Extended results for: The balance of nature relies on response diversity for stability

Til Hämmig, Francesco Polazzo

13 November, 2024